COMMENT ⊗   VALID 00002 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	\magnification\magstephalf
C00015 ENDMK
C⊗;
\magnification\magstephalf
\input macro.tex
\def\today{\ifcase\month\or
  January\or February\or March\or April\or May\or June\or
  July\or August\or September\or October\or November\or December\fi
  \space\number\day, \number\year}
\baselineskip 14pt
\rm
\line{\sevenrm a66.tex[106,phy] \today\hfill}

\bigskip
\line{\bf Monte Carlo Methods.\hfill}

There is no part of a computer that resembles dice, cards, a spinning coin,
or any other such device prized for its randomness.  Still, at times we
want to use a computer to study the behavior of systems that do contain
random elements.  The telephone company tries to ensure that its circuits
will accomodate not only average traffic, but most of the occasional
``random'' peaks.  Similar studies can be made of vehicle traffic, customers
queueing (yes, it's an English word) in a bank, or evolution of bird
species.  Often there seems to be no better way to analyze such partially
random systems than to simulate them: to write computer programs that
behave as they do, including their random elements.

Just as all polynomials can be built up by addition and multiplication, so
can most random phenomena be built up from a process that picks a random
number.  Ordinarily, the random number is picked from a certain range, its
possible values in that range are equally spaced, and they are all equally
likely.  A process that does this we call a {\sl random number generator\/}.  It
might pick an integer between~0 and {\tt MAXINT}, or a real number between~0 
and~1; either is useful.  In a computer, a random number generator will
typically be a function or procedure that returns a value according to the
rules above.

Let's assume we have a random number generator procedure {\tt RANDOMREAL}, which
gives its variable argument a real number value between~0 and~1.  The
chance that the value lies between {\tt A} and~{\tt B}, 
if $0≤{\tt A}<{\tt B}<1$, is roughly ${\tt B}-{\tt A}$.
Using  {\tt RANDOMREAL},
we can simulate a dice roll. After the procedure call
{\tt RANDOMREAL(X)}, we can say the number on the die is~1 if $0≤{\tt x}<1/6$, 
2~if
$1/6≤{\tt X}<2/6$, etc.  
The formula {\tt TRUNC(X/6)+1} gives the number on the die; it
is equally likely to have any value from~1 to~6.

A more complicated random system to simulate is an atom of a radioactive
element.  Say it has a 1\%\ chance of decaying in each second, if it has not
already done so.  That is, the chance it decays in the first second is
0.01.  The chance it decays in the second second is 0.01 of the remaining
99/100 of the time, or $0.01\times  0.99$.  The chance it decays in the 
$n$th~second is then $0.99↑{n-1}\times  0.01$.  
If our random number generator
{\tt RANDOMREAL(X)} gives us a real number {\tt X} between~0 and~1, 
we can compute a time $t$ for the atom to decay by setting 
$$\vcenter{\halign{$\rt{#}$\qquad&\lft{#}\cr
t=1&if $0.99≤X<1$;\cr
t=2&if $0.99↑2≤X<0.99$;\cr
t=3&if $0.99↑3≤X<0.99↑2$, etc.\cr}}$$

Taking logarithms, we see
$$\vcenter{\halign{$\rt{#}$\qquad&\lft{#}\cr
t=1&if $1≥\log X/\log 0.99 > 0$\cr
t=2&if $2≥\log X/\log 0.99 > 1$\cr
t=3&if $3≥\log X/\log 0.99 > 2$, etc.;\cr}}$$
we can use the formula t= truncate $(\log X/\log 0.99) + 1$.

To simulate an event that has probability {\tt P} of occurring, generate a random
real number~{\tt X} in the range {\tt 0..1}; 
the event occurs if {\tt X<P}.  Alternatively,
generate a random integer {\tt X} in the range {\tt 0..MAXINT}; the event occurs if
{\tt X < P * MAXINT}.  

There is no Pascal standard random number generator.  Many Pascal systems
include one and users should ordinarily prefer such a ``certified'' random
number generator to a homemade one.  The typical random number generator
requires a starting value, called a {\sl seed\/}, which determines the sequence of
subsequent values.  If the same program is run a second time with the same
seed, the same results will be found, so to carry out statistically
independent executions, each time the program is run, a new seed should be
used.  The seed should be recorded, so that if bugs are found or suspected
the program can be run again under identical conditions.

A typical programming cliche for using a seeded random number generator is:

\smallskip\halign{\qq\qq\lft{\tt #}\cr
	READ(SEED);\cr
	WRITELN('SEED=', SEED);\cr
	X:=SEED;\cr
	WHILE ... DO\cr
\qq		BEGIN\cr
\qq		RANDOM(X);\cr
\qq		Process random number X\cr
\qq		END\cr}

\smallskip
It is not ordinarily  possible to generate truly random numbers on a
correctly operating computer; we satisfy ourselves by generating sequences
that meet most statistical tests for randomness, like the requirement that
in a large trial, nearly equal fractions of the random numbers fall in each
percent of the possible range.  The typical method, using an integer~{\tt X}, is

\smallskip\halign{\qq\qq\lft{\tt #}\cr
	X:= (X * A)MOD M\cr}

\smallskip\noindent
where {\tt A} and {\tt M} are carefully selected constants.  A~built-in 
random number
generator usually uses double precision integer arithmetic, so it can not
easily be programmed in standard Pascal.  To demonstrate the pattern of
behavior, here is the sequence of random numbers if ${\tt A}=23$, 
${\tt M}=100$, and the seed is~47:

\smallskip\halign{\qq\qq\lft{\tt #}\cr
81 63 49 27 21 83 09 07 61 03\cr
69 87 01 23 29 67 41 43 89 47\cr
81 63 {\rm etc.}\cr}

\smallskip
Among the flaws of this choice, the sequence repeats after generating only
twenty numbers, and it generates only numbers that end in 1,3,7, or~9.  If
the seed had been~25, the sequence would have been 25 75 25 75 25, etc.  By
choosing ${\tt A}=??$, ${\tt M}=??$, however, we could guarantee that every number
between~1 and~??  would be generated  once before repeating, for any seed
in that range.

If you must build your own random number generator, see Knuth, [fill in reference]
for some rules of thumb which ensure good performance.

Here is a Pascal program using an integer random number procedure to
simulate a dice game.  Player~1 rolls a die once, then player~2 rolls it
twice, player~1 rolls it three times, etc., until someone rolls a~6.  The
program simulates the game a number of times to estimate the probability
that {\tt A} wins the game.

\smallskip\halign{\qq\qq\lft{\tt #}\cr
     CONSTANT C=MAXINT/6\cr
     READ(SEED, NUMTRIALS)\cr
     X:=SEED;\cr
     WINS:=0;\cr
     FOR I:=1 TO NUMTRIALS DO\cr
\qq	BEGIN\cr
\qq	PLAYER:=1;\cr
\qq	CHANGE:=1;\cr
\qq	DIE:=0;\cr
\qq	COUNTER:=0;\cr
\qq	WHILE DIE <> 6 DO\cr
\qq\qq	     BEGIN\cr
\qq\qq	     RANDOM(X);\cr
\qq\qq	     DIE:=X DIV C + 1;\cr
\qq\qq	     COUNTER:= COUNTER + 1;\cr
\qq\qq	     IF COUNTER = CHANGE THEN\cr
\qq\qq\qq	BEGIN\cr
\qq\qq\qq	COUNTER:=0;\cr
\qq\qq\qq	CHANGE:=CHANGE + 1;\cr
\qq\qq\qq	PLAYER:=3 - PLAYER;\cr
\qq\qq\qq	END\cr
\qq\qq	     END;\cr
\qq	IF PLAYER =1 THEN WINS:= WINS + 1\cr
\qq	END;\cr
	WRITELN(`PLAYER 1 WINS', WINS, `OUT OF', NUMTRIALS);\cr
	WRITELN(`FOR AN ESTIMATED PROBABILITY', WINS/NUMTRIALS);\cr
	WRITELN(`SEED WAS', SEED)\cr}

\bigskip
\line{\bf Exercise:\hfill}

Simulate an island populated by wolves and rabbits.  The island is square,
and looks very much like a 50-by-50 array of regions.  Each region of the
island may be empty, or contain a wolf or a rabbit.  At each day, each wolf
has a 25\%\ probability of moving to any adjacent wolf-free region, eating
any rabbit found there.  Rabbits have a 25\%\ probability of moving to any
adjacent empty region.  Rabbits have a 5\%\ probability per day of
reproducing [FILL IN].  A~wolf starves if it has not eaten a rabbit
for ??  days.  Wolves have a   ??  probability per day of reproducing.

Simulate the population of the island for 1000 days, printing a map showing
locations of the animals every $200↑{\rm th}$ day, and tabulating the populations at
10~day intervals.  Initially, place about 10~wolves and 500~rabbits at
random. 








\bigskip
\parindent0pt
\copyright 1984 Robert W. Floyd

First draft October 31, 1984

\bye